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Abstract 

Rate coefficients of the Arrhenius-Neel form are calculated for thermally activated magnetic 
moment reversal for dual layer exchange-coupled composite (ECC) media based on the Langer for¬ 
malism and are applied to study the sweep rate dependence of MH hysteresis loops as a function 
of the exchange coupling / between the layers. The individual grains are modelled as two exchange 
coupled Stoner-Wohlfarth particles from which the minimum energy paths connecting the mini¬ 
mum energy states are calculated using a variant of the string method and the energy barriers 
and attempt frequencies calculated as a function of the applied field. The resultant rate equations 
describing the evolution of an ensemble of non-interacting ECC grains are then integrated numer¬ 
ically in an applied field with constant sweep rate R = —dH/dt and the magnetization calculated 
as a function of the applied field H. MH hysteresis loops are presented for a range of values / for 
sweep rates 10 5 Oe/s < R < 10 10 Oe/s and a figure of merit (FOM) that quantifies the advantages of 
ECC media is proposed. MH hysteresis loops are also calculated based on the stochastic Landau- 
Lifshitz-Gilbert equations for 10 8 Oe/s < R < 10 10 Oe/s and are shown to be in good agreement 
with those obtained from the direct integration of rate equations. The results are also used to 
examine the accuracy of certain approximate models that reduce the complexity associated with 
the Langer based formalism and which provide some useful insight into the reversal process and 
its dependence on the coupling strength and sweep rate. Of particular interest is the clustering of 
minimum energy states that are separated by relatively low energy barriers into “metastates”. It is 
shown that while approximating the reversal process in terms of “metastates” results in little loss 
of accuracy, it can reduce the run time of a Kinetic Monte Carlo (KMC) simulation of the magnetic 
decay of an ensemble of dual layer ECC media by 2 ~ 3 orders of magnitude. The essentially exact 
results presented in this work for two coupled grains are analogous to the Stoner-Wohlfarth model 
of a single grain and serve as an important precursor to KMC based simulation studies on systems 
of interacting dual layer ECC media. 
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I. INTRODUCTION 


After many decades, the Landau-Lifshitz-Gilbert (LLG) equation continues to provide 
the foundation for micromagnetic modeling of the dynamic evolution of granular magnetic 
material, with an increasing number of applications devoted to the study of the effects of ther¬ 
mal fluctuations. 12 The LLG equation is commonly used to study granular recording media 
but it is limited to the study of phenomena over relatively short time scales (ms). Modern 
exchange-coupled composite (ECC) recording media is composed of a high anisotropy ‘hard’ 
layer exchange coupled to one or more lower anisotropy ‘soft’ layers.^ 9 One of the most im¬ 
portant applications of micromagnetic modeling is the characterization of recording media 
through MH hysteresis loops. Often, model parameters are determined by fitting results to 
experimental data. Inconveniently, experimental hysteresis loops typically require minutes 
to hours to complete, time scales that are outside the range of standard LLG simulations. 
In addition, the thermally activated decay of recorded bits requires a micromagnetic model 
that is valid over much longer time scales (years). Various scaling arguments, based on the 
Arrhenius-Neel law, have been proposed as a means to extrapolate LLG results to longer 
times which appear useful for older single layer type recording media.® Thermally activated 
processes in ECC media, however, are more complex and the simple scaling arguments 
appear to break down in this case.®! 

For the purpose of studying long-time scale micromagnetics governed by thermally acti¬ 
vated processes, Kinetic Monte Carlo (KMC) methods have proven useful. 12 — In Ref. [T71 
a KMC algorithm to study long-time scale thermally activated grain reversal of single layer 
recording media was described. The Arrhenius-Neel expression for the rate coefficients be¬ 
tween the minimum energy states of the individual grains was used to calculate the time 
between successive reversals. The minimum energy states and the energy barriers sepa¬ 
rating them were calculated using a modified version of the Wood analytic expression for 
single Stoner-Wohlfarth particle^ 1 ® (SWPs) which includes the effective exchange and mag¬ 
netostatic fields from neighbouring grains. For weakly interacting recording media, the 
effective field approximation appears valid. For the attempt frequency, the temperature 
and field dependent formula of Wang and Bertram, 19 based on a single energy barrier was 
used. This algorithm was subsequently used to study the magnetic MH hysteresis loops of 
high anisotropy magnetic recording media at both short and long time scales over a wide 


3 


range of temperatures relevant to heat assisted magnetic recording. 20 Good agreement be¬ 
tween the KMC results and those from LLG simulations at relatively short time scales was 
demonstrated. 

In Ref. [2U this KMC algorithm was applied to study MH hysteresis loops of dual layer 
ECC recording media at finite temperature and long time scales in which the effect of 
ECC interlayer exchange coupling was treated in the same way as the intralayer exchange 
coupling by means of an effective field. ECC media for practical applications has a relatively 
strong exchange coupling and it is not evident that treating the intralayer coupling through 
an effective field is a good approximation 8 for this purpose as it ignores the correlated 
nature of the rotation of the layers in the reversal process. In addition, the expression used 
for the attempt frequency is based on a single energy barrier, is unlikely to be valid for 
multi-layer ECC media at relevant (moderate) coupling strengths. The absence of simple 
Arrhenius-Neel-type scaling between thermal and temporal effects for ECC media supports 
these conclusions. 11 

In the present work, we study the reversal process for two interacting magnetized grains 
which treats the correlated reversal process in a more systematic way. The approach includes 
the complete set of minimum energy states for the ECC grains, while the calculation of the 
rate coefficients, based on the Langer formalism,takes into account the complex minimum 
energy paths (MEPs) connecting them. The resultant energy barriers and attempt frequen¬ 
cies provide a comprehensive treatment of reversal process that is applicable to both weak 
and strong interlayer coupling. 

In order to study statistical effects, we also consider an ensemble of non-interacting 
exchange coupled dual layer grains. This has the benefit that we can describe the evolution 
of the ensemble from some initial distribution of states in terms of a system of rate equations 
that can be integrated numerically. In particular we present a series of MH hysteresis loops 
calculated at constant sweep rate over a range of coupling constants to examine the effect of 
the exchange parameter and sweep rate on the MH hysteresis loops. This work compliments 
and extends previous studies of dual-grain-reversal energy landscapes® and formulations of 
the dual-grain attempt frequency,^ 24 and serves as an precursor to our formulation and 
application of the combined MEP-KMC algorithm to study interacting NxNx2 ECC thin 
films which includes both magnetostatic and intralayer exchange interactions. 25 

In the following section, we discuss the energy landscapes for a dual layer grain, which 
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we model as a system consisting of two coupled Stoner-Wohlfarth particles (SWPs) in a 
magnetic field. We briefly outline how the rate coefficients may be calculated based on 
the Langer formalism for such a system of exchange coupled SWPs in the strong and weak 
coupling regimes. The details of the Langer formalism are presented in the Appendix. In 
Sec. |III[ we show how the equations may be integrated and the rate coefficients determined 
from the MEP calculated using a variant of the so-called “string method” E5HZS an d a series 
of Addd hysteresis loops are presented for various sweep rates and couplings. A figure of 
merit to assist in the evaluation of the benefits of ECC media is proposed based on the 
ratio of the switching field and the energy barriers and is calculated as a function of the 


exchange coupling for several different sweep rates. In Sec. |IV[ we compare the Addd hys¬ 
teresis loops obtained from the rate equation approach with those obtained using stochastic 
LLG simulation over a range of sweep rates, where both approaches should be valid. The 
good agreement between the results from the two methods gives confidence that the results 
obtained from the rate equation approach, which can be extended to very low sweep rates, 
are essentially of the same quality as those obtained from stochastic LLG, which is restricted 
to very high sweep rates. 


In addition to the MEP based calculations of the thermally activated dual-grain reversal, 
we also examine the validity of two important approximation schemes. The first, based on 
approximations to the MEP allow us to obtain analytical expressions for the energy barrier 
and attempt frequency in the strong and weak coupling limits, respectively. A comparison 
of results obtained based on this scheme with MEP calculations are presented in Sec. |V| and 
show good agreement over a range of couplings and sweep rates. The second approximation 
method exploits the fact that in dual layer ECC media, the rate coefficients calculated from 
the MEPs separating pairs of energy minima can differ by orders of magnitude. A direct 
consequence of this is that pairs of minima will equilibrate on time scales that are significantly 
shorter than the time taken to complete a single sweep, or portion of a sweep. In such cases 
it is possible to combine the two minimum energy states into a single “metastate”. In 
Sec. |VTJ we show how this approach can be used to reduce a 4-state model to an equivalent 
2-metastate representation that gives essentially the same results. While the difference 
in computation time required to solve the rate equations for the 4-state model and the 
equivalent 2-metastate model is negligible, the same is not true when we use KMC to 
calculate Addd hysteresis loops that include the interlayer interaction between the layers. 
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In this case, the large variation in rate coefficients causes the KMC algorithm to slow to 
a crawl, giving rise to what we refer to as “stagnation”. In Sec. VII we illustrate the 
effects of stagnation by applying the KMC method to compute magnetization decay for 
an ensemble of non-interacting ECC grains. It is shown that by combining certain pairs 
of states into metastates, the calculation speeds up by a factor of 600 with negligible loss 
in accuracy. KMC studies on dual layer ECC grains that include the magnetostatic and 
intralayer exchange interactions show that the clustering of minimum energy states separated 
by low energy barriers into metastates to avoid the effects of stagnation is critical to the 
successful application of KMC to systems of interest in magnetic recording media. 25 


II. ENERGY LANDSCAPES AND RATE EQUATIONS IN ECC MEDIA 

ECC media consists of magnetic grains with different layers of varying anisotropy strength 
and moderate exchange interactions between these layers.® 3 8 The desired effect is to be able 
to use the very strong anisotropy of a hard layer to enhance the grains thermal stability 
and hence prevent data loss due to thermally activated grain reversal. The hard layer is 
then exchange coupled to a layer with lower anisotropy, the soft layer. The soft layer will 
respond more readily to a switching field and the exchange coupling interaction will make 
switching the hard layer easier. The result is a thermally stable grain that can be reoriented 
by using an applied field at lower magnitudes than would be required if just the hard layer 
were presented. 

In this work, we consider an ensemble of two exchange coupled grains. In order to focus 
on the role of the exchange coupling between the layers and to allow a semi-analytical 
treatment of the system, we neglect the lateral exchange interaction between the grains and 
the magnetostatic interaction. 

The grains are cubic with a side length a=6 nm stacked along the z-axis. The dimensions 
of the grains are such that each may be treated as a single domain ferromagnet and may be 
modelled as two exchange coupled SWPs which we label as a and b. The energy of a single 
grain is therefore written in terms of the normalized magnetization vectors rhi = Mi/Mi, 

E = - K a v a (m a ■ h a ) 2 - K b v b ( m b ■ h b ) 2 - IA(m a ■ m b ) 

- id 0 H ■ ( M a v a m a + M b v b rh b ), (1) 
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where H denotes the applied held and A7, Vi, hi and M* denote the anisotropy con¬ 
stant, volume, anisotropy axis and the saturation magnetization of the i th grain, respec¬ 
tively. The grains we consider are comprised of a soft layer with M a =4.0xl0 5 A/m and 
K a — 1.5xl0 5 J/m 3 , and a hard layer with M b =5Ax 10 5 A/m and A/,=3.0xl0 5 J/m 3 . The 
layers are coupled through ferromagnetic exchange expressed in terms of the coupling con¬ 
stant I and interfacial area A = a 2 . Assuming that both the anisotropy axis h* and the held 
H are aligned perpendicular to the plane, then the SWP energy for a single grain may be 
written in spherical coordinates as, 

E = - K a v a sin 2 6 a - K b v b sin 2 6 b - fi 0 H (M a v a sin 6 a + M b v b sin 6 b ) 

— IA (sin 9 a sin 9 b + cos (4> a — <$> b ) cos 9 a cos 9 b ) , (2) 

where 6 a and 6 b denote the polar angles measured relative to the xy plane and (j) a and (j) b 
denote the azimuthal angles associated with the grains a and b , respectively. 

To understand how the energy of a grain depends on the variables {6 a , 4> a , 6 b , (ft b }, we note 
hrst that it is invariant under rotation about the z-axis and thus depends only on the three 
independent variables {9 a , 6 b , (j) b — cj) a }. Also since the exchange coupling between the layers 
is positive, the energy is minimized when (p a = (p b , we therefore find it useful to plot the two 
dimensional subspace defined by 0 a = (j) bl which refer to as the “minimum energy surface”. 
We consider four specific cases in some detail corresponding to I = 2.0 x 10~ 3 J/m 2 and 
/ = 0.5 x 10 -3 J/m 2 for both H = 0 and y 0 H = 4 kOe. 


A. Strong Exchange Coupling 


In Fig. 1(a) the contour plot of the minimum energy surface over the range —7 t/2 < 9 a < 
7t/ 2 and —7 t/2 < 9 b < 7t/ 2 for H = 0 and I = 2.0 x 10~ 3 J/m 2 is presented. The energy 
landscape shows two minima corresponding to two stable states with the magnetic spins 
aligned ferromagnetically along the z-axis. We refer to this as the strong exchange coupling 
regime and denote the two minimum energy states {9 a , 9 b } = {— 7t/2, —7t/ 2} and {7 t/2, 7t/ 2} 
as <7i and <r 4 , respectively. We note that in the absence of a held, the energy landscape 
is symmetric under spin inversion and hence the minimum energy states are degenerate 
Ei = E 4 . 


The contour plot of the minimum energy surface of a grain is presented in Fig. 1(b) for 
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H — 4 kOe. The minimum energy landscape again shows two minima located at cq and <74 
however, because of the applied field, the energies are no longer degenerate and E A < E\ so 
the system now has one stable minimum at (74 and a metastable minimum at eri. 



FIG. 1. The contour plot of the minimum energy surface for the strong coupling case 
7=2.0 xlO 3 J/m 2 for (a) 77=0 and (b) 77=4 kOe. The black lines indicate the boundary sep¬ 
arating the two basins of attraction fh and O 4 which we denote by The red lines indicate 
the MEP connecting the minimum energy states o\ and < 74 . The lines cross at the saddle point 
indicated by S 14 . 


Associated with each of the local minimum energy states cr a is a basin of attraction 
defined as the region of phase space comprising the states that evolve asymptotically to the 


state a a . We denote the basin of attraction associated with the state a a as tt a . Figure 1(a) 


and Fig. l(b)| show the boundary separating the two basins of attractions fli and Q 4 which 
we denote by T 14 . 

The probability distribution of the grains in phase space is given by the probability 
density p(x,t), where x denotes a vector that specifies a point in phase space in terms 
of some generalized coordinates (e.g. (x 1 , x 2 , x 3 , x 4 ) = (9 a ,(pa,9b,(pb))- The evolution of 
the probability density is given by the Fokker-Planck equation (FPE). For the energy and 
time scales of interest to us, the system will be in local equilibrium. Local equilibrium 
assumes that, except for a narrow crossover region A 14 that runs along the boundary r 14 , 































the probability density p(x, t) is given by the Boltzmann distribution, 

p(x,t) c a (t) exp forx€fi a -Ai 4 , (3) 

where (except in the case of thermal equilibrium) C\^ c^- The probability p a that a grain 
in the ensemble is located in is therefore given by, 

Pa{t) = C a (t ) f exp [- \ dtt = c a (t)Z a . (4) 

In the crossover region A 14 , the system is not in equilibrium and the probability density 
is given by the more general form p(x,t ) = c(x,t) exp(—E(x,t)/kBT), where the crossover 
function c(x,t) is obtained from the FPE and interpolates between the coefficients C\ (f) 
and 04 (f) defined by Eq. ([3]). The inhomogeneous nature of c(x,t) in the crossover region 
gives rise to a net flux of probability across the boundary that is driven by the thermal 
fluctuations. For the energy scales of interest, this probability flux is concentrated at the 
point of minimum energy on the boundary T^. This point, which we denote by S 14 , is a 
saddle point with dE/dx^ = 0 and a Hessian matrix \\d 2 E/dx^dx v \ \ that has two positive 
eigenvalues, one negative eigenvalue and a zero eigenvalue (the latter arising as a consequence 
of the rotational symmetry of the energy about the axis perpendicular to the plane). 

The rate at which the particles in the ensemble make the transition from cr a — > <jp may 
be expressed in terms of the rate constant as, 


>7? ^a/3Po(^); 


( 5 ) 


where the rate constants r a p are of the Arrhenius-Neel form, 


r a p = f a p exp 


f _ AE a p\ 


( 6 ) 


where the energy barrier AF a ^ = E(s 14 ) — E(cr a ) and the attempt frequency, / Qj g, may be 
calculated from the crossover function C(X,t ) in the neighbourhood of the saddle point s a p, 
and may be expressed as, 


_ « 0 j g(s) 7 B r , , / 1 V 1 V 2 WU m 

1 + al\Jg(a)m U K y 2nk B T |AiA 2 A 4 | ’ 

where «o is the damping constant, 7 b is the gyromagnetic ratio, m = M a v a + Mi,Vb, X t 
and rji are the eigenvalues of the Hessian matrix \\d 2 E/dx^dx v \\ calculated at the saddle 
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point, s a p, and the minimum energy state, a a , respectively, k -1 / 2 characterizes the width 
of the crossover region A a p in the vicinity of the saddle point s a p, and the quantities g(s), 
g(a), and Gu are related to metric associated with the particular coordinate system (or 
systems) used in the derivation. The details of the derivation of Eq. ([7]) are presented in the 
Appendix. 

The time dependence of the probabilities p a (t) can be calculated from the rate equations, 


dpi 

dt 

dp 4 

dt 


-tiaPi + r 41P4 
-r 41 p 4 + r 14 pi. 


( 8 ) 

(9) 


B. Weak Exchange Coupling 


Figures [2(a)| and 2(b) show the corresponding energy landscapes for the case / = 0.5 x 
10 ~ 3 J/m 2 for H = 0 and 4 kOe. Both cases have four minimum energy configurations 
corresponding to the ferromagnetically aligned states {9 a ,9 b } = {q=7r/2, =F 7 t/ 2 } and the an- 
tiferromgnetically aligned states {9 a ,9 b } = {± 7 t/ 2 , =F 7 t/ 2 }. We denote the minimum energy 
states {=F7 t/ 2, q= 7 r/ 2 } by a 4 and < 74 , as in the strong coupling case discussed above, and 
the antiferromagnetic states { 7 t/ 2 , — 7 t/ 2 } and {— 7 t/ 2 , 7 t/ 2 } as 02 and cr 3 , respectively. We 
refer to this as the weak exchange coupling regime. As before, in the absence of an applied 
held, the system is invariant under a spin inversion and we have the following degeneracies 


Ei = E 4 and E 2 = E 3 . Figures 2(a) and 2(b) also show the basin boundaries T a p separating 
the basins of attraction associated with the energy minima cr a . 

For the energy and time scales of interest, the system will be in local equilibrium and we 
can therefore define the probabilities that a grain is located in Q a as, 


Pa(t) = c a (t ) f exp d£l = c a (t)Z a (t ), (10) 

and the probability flux T a ^.p between the basins of attractions will be concentrated 

at the saddle point s a p and may be expressed in terms of the rate constants r a p of the forms 
given by Eq. ([5]). The formalism presented in the Appendix applies equally well to the 
systems with multiple energy minima. The rate equations given by Eqs. (J8]) and ([9]) for 
strong coupling regime may then be written to include the case of multiple (i.e. more than 
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FIG. 2. The contour plot of the minimum energy surface for the weak coupling case I = 0.5 x 
ICG 3 J/nr 2 for (a) H = 0 and (b) H = 4 kOe. The black lines indicate the boundary separating 
each pair of basins of attractions {fia, and is denoted by T a p. The red lines indicate the MEP 
connecting the four minima. The lines cross at the saddle points s a p. 


two) minima as, 


dp a 

dt 


N s 

Y ( r c*pPa - rga'Pd). 

13=1 


( 11 ) 


In applying the above formula we note that r aa = 0 and that the number of minimum energy 
states, N s , will depend on the strength of the applied field, ranging from 1 to 2 in the strong 
coupling regime and from 1 to 4 in the weak coupling regime. 


III. THE MINIMUM ENERGY PATH AND THE EVALUATION OF MH HYS¬ 
TERESIS LOOPS 

To evaluate the MH hysteresis loop for a layer of non-interacting ECC grains using the 
rate equations given by Eq. 0. we consider that at some initial time, t = tj, the system is 
fully saturated p\{t = U) = 1 in a large positive applied field with only one minimum energy 
state <7i. The held is then reduced at a constant rate dH/dt = —R until the system is again 
fully saturated in the opposite direction at time tf, p±{t = tf ) = 1. Since the rate of change 
of the applied held is constant, we have that dp a /dt = —Rdp a /dH and the rate equations 
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may be written as, 


dp a (H) 

dH 


N s 


RJ2 MH)p a (H) - r Ba (H)p B (H)). 


( 12 ) 


3=1 


Integrating these equations with the initial condition p±(H = Hi) = 1 yields the probabilities 
p a (H) from which we can then compute the magnetization as a function of H, 


m(H) = m aPa(H), 


(13) 


where m a denotes magnetic moment of a grain in state cr a . 

Calculating the rate constant r a p as a function of H requires that we determine the 
location of the minimum energy states and the saddle point located on the boundaries 
separating their basins of attraction for each held value. For this simple example, locating 
the energy minima is straightforward. How best to determine the location of the saddle 
points is less obvious. One technique is to compute the MEP that connects the two minima 
a a ap. This may be done numerically by discretizing an initial guess of the MEP and 
allowing the points to relax until the derivatives of the energy perpendicular to the tangent 
line at each point of the path are zero. Two methods that successfully implement this 
scheme are the Nudged Elastic Band (NEB) method and the string method . 26 28 In the 
present work, we have used a variant of the string method to calculate the MEPs. MEPs for 
both the strong coupling (/ = 2.0 x 10 -3 J/m 2 ) and the weak coupling (/ = 0.5 x 10 -3 J/m 2 ) 
regimes are shown in Figs, [l] and [2] for both H =0 and H =4 kOe. It should be noted that 
not every pair of energy minima are directly connected by an MEP. For such cases r a p = 0. 


Parametric plots of the energy E(6 a ,6b ) along the MEP are shown in Figs. 3(a) and 3(b) 


for / = 2.0 x 10 -3 J/m 2 and in Figs. 4(a) - 4(d) for / = 0.5 x 10 ” 3 J/m 2 . The figures also 
show parametric plots for the energy along the initial path. For the strong coupling regime 
(/ = 2.0 x 10 ~ 3 J/m 2 ), the initial path used is given by (0,9), while in the weak coupling 
regime (I = 0.5 x 10 ~ 3 J/m 2 ), where there are up to four MEPs, the initial paths used were 
(6,— 7 t/2 ), (ti/2,0), (—ti/2,0), and (0,n/2) where — n/2 < 0 < tt/2. The saddle point is 
located at the point of the peak energy on the MEP. 


The integration of Eq. (12) proceeds as follows, the minimum energy states for several 


values of H over the range —5 kOe < poH < 2 kOe and the MEPs joining them are 
determined. The saddle points are located at the point of maximum energy on the MEP. 
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FIG. 3. A plot of the energy along the length of the parametric MEP (red line) and the initial 
path used to generate it (black line) for 7=2.0xlCF 3 J/m 2 (a) 77=0 and (b) 77=4 kOe. 


Once the minimum energy states and the saddle points have been determined, the non¬ 
zero rate constants r a/ g(77) are calculated using the expression given by Eqs. @ and 0 
at these selected values of 77. The rate coefficients for values of 77 intermediate between 


these discrete values are then determined by interpolation. The rate equation (12) is then 
be solved numerically using Mathematica. 

The range of sweep rates chosen corresponds approximately to time scales involved in 
experimental MH hysteresis loop measurements, R ~ 10 3 Oe/s, to magnetic recording rates, 
R ~ 10 10 Oe/sP Figure [5] shows the calculated MH hysteresis loops at T=300K with 
a 0 = 0.1 for different exchange coupling values 7 = 2.0,1.5,1.0, 0.5, 0.25 and 0.1 x 10 _ 3 J/m 2 , 
respectively. MH hysteresis loops are calculated at all the different sweep rates, but only 
the range 10 5 Oe/s < R < 10 10 Oe/s are shown in these figures. 

From these results, the expected trend of the coercivity H c decreasing at slower sweep 
rates can be observed. In addition, there is little difference between the strong coupling 
cases of 7=2.0xl0 -3 J/m 2 and 7=1.5x 10 ~ 3 J/ 111 2 . Moderate coupling 7=1.Ox 10 ~ 3 J/m 2 
and 7=0.5x 10 -3 J/m 2 represents a crossover regime between the two grains acting as a 
single grain, and the two grains responding quasi-independently. Here the hysteresis loops 
are quite sensitive to the coupling 7. Weak coupling is clear in the case of 7 = 0.1xl0 ~ 3 J/m 2 
at the fast sweep rate, where the plateau indicates that the soft grain switches Erst. 

From the hysteresis loops, we can extract the nucleation held H n = H(M/M s = 0.95), the 
coercive held T7 C = H(M/M S = 0.0), and the saturation held H s = H(M/M S = —0.95 ). 20 
Figure [ 6 ] shows these extracted values of 77„, 77 c , and H s as a function of 7 for different 
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FIG. 4. A plot of the energy along the length of the four parametric MEPs (red and blue lines) 
and the initial paths used to generate them (black and green lines) for 1= 0.5xl0 -3 J/m 2 (a) H =0 
and (b) H =4 kOe. 


sweep rates. Although the nucleation field exhibits monotonic decrease with increasing / 
and R , both H c and H s show clear minima at weak to moderate coupling values in the cases 
of the faster sweep rates. 


IV. COMPARISON WITH LLG 


As mentioned above and discussed in more detail in the Appendix, the calculation of 
the rate coefficients follows from the FPE, which can be derived from the stochastic LLG 
equation . 22 29 30 In fact the derivation of the FPE from the stochastic LLG equation imposes 
non-trivial requirements on the integration schemes that can be used to solve the stochastic 
LLG equation. It is therefore interesting to compare the MH hysteresis loops obtained from 


stochastic LLG and those obtained in Sec. Ill Previous comparisons for interacting grains 
where the rate equations have been solved using both Kinetic Monte Carlo (KMC ) 17 20 and 
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FIG. 5. MH hysteresis loops calculated at T=300 K by direct integration of the rate equations for 
different sweep rates R: (a) /=2.0xl0~ 3 J/m 2 , (b) /=1.5xl0 -3 J/m 2 , (c) /=1.0xl0 -3 J/m 2 , (d) 
1= 0.5xl0 -3 J/m 2 , (e) /=0.25xl0 -3 J/m 2 , and (f) 1= O.lxlO -3 J/m 2 . 


stochastic LLG have shown good agreement between the two approaches over a limited range 
of sweep rates (10 8 Oe/s < R < 10 10 Oe/s). The range of sweep rates over which we might 
expect good agreement between the two approaches is limited by the fact that LLG results 
are only accessible within a reasonable amount of simulation time for R > 10 8 Oe/s while 
for R > 10 10 Oe/s, the Arrhcnius-Neel expression for the rate coefficient, that serves as a 
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R = 10 x (Oe/s) 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 


I(10' 3 J/m 2 ) 


R = 10* (Oe/s) 



FIG. 6. (a) The nucleation field H n , (b) coercivity H c , and (c) saturation field H s extracted from 
Fig. H as a function of I at different sweep rates. 


basis for the KMC algorithm, breaks down, as it does not fully capture the spin dynamics 
of the reversal process. 


MH hysteresis loops obtained from LLG simulations for a system of 16 x 16 non¬ 
interacting, exchange coupled dual layer grains using the same parameters detailed in Sec. [TT] 
are presented in Figs. 7(a)| - [7(cj| together with loops obtained by the MEP method. The time 
step used was 2 ps and the integration was performed using the Runge-Kutta fourth order 
method based on a quaternion representation of the rotations with the damping parameter 


set at cto = 0.1. The simulations performed at T = 300K. The comparison shown in Fig. 7(a) 


for 7=2.0x10 3 J/m 2 , Fig. 7(b) for 7=0.5x10 3 J/m 2 , and Fig. 7(c) for 7=0.1x10 3 J/m 2 
indicates a very good agreement between the two methods at all sweep rates. 
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H (kOe) 

FIG. 7. Comparison of MH hysteresis loops from the MEP method (solid curves) and stochastic 
LLG (dashed curves) at different sweep rates for (a) 7=2.0xl0 -3 J/nr 2 , (b) 7=0.5 xl0 ~ 3 J/m 2 , and 
(c) 7=0.1xl0 ~ 3 J/m 2 . 

V. FIGURE OF MERIT FOR ECC MEDIA 


The benefit of coupling hard and soft layers can be quantified in a figure of merit (FOM), 
which is the ratio of a measure of the thermal stability and the field required to switch 
the grain magnetization . 3 8 This can be defined as the ratio between the energy barrier Eb 
(thermal stability) and saturation field (switching energy) at a particular sweep rate as, 


rp 

FOM =---. (14) 

fx 0 H s (M a v a + M b v b ) 

For strong coupling, Eb is given by the zero field energy barrier between the minimum 
energy of state cri and the saddle point along the path to the minimum energy of state 
a. 4 , while for weak coupling, Eb is the zero field energy barrier between the minimum en¬ 
ergy of state < 71 and the saddle point along the path to the minimum energy of state cr 2 . 


A larger FOM is the goal for ECC-type media. The results shown in Fig. 6 (c) indicate 


17 



































that increasing /, for small /, will decrease the saturation held which makes switching the 
magnetic moment easier. On the other hand, increasing / will increase the energy bar¬ 
rier which enhances the thermal stability (not shown). Fig. [ 8 ] shows the FOM at three 
sweep rates (R = 10 6 , 10 8 , and 10 10 Oe/s), and the optimal value of / can be easily ob¬ 
tained from the graph: J op (10 6 Oe/s)~0.2x 1CT 3 J/m 2 , J op (10 8 Oe/s)~0.35x 10 ~ 3 J/m 2 , and 
J op (10 10 Oe/s)~0.50x 10 “ 3 J/m 2 . These results suggest that weak to moderate coupling is 
preferred and that there is a strong dependence on sweep rate. Large FOM values at smaller 
sweep rates may not be realized at larger sweep rates, and optimal coupling strengths esti¬ 
mated on the basis of experimental MH hysteresis loops obtained at slow sweep rates may 
thus not be the optimal value at recording time scales. 



FIG. 8. Figure of merit calculated by Eq. (14) for three sweep rates (10 6 , 10 8 , and 10 10 Oe/s) 


VI. APPROXIMATION SCHEMES 


In this section we describe some approximation schemes which allow simplification of the 


rate equations used in Sec. Ill, not only making calculations less onerous but also, in certain 
cases, allowing for analytical solutions. Comparisons with the exact rate equations show 
that for certain regions of parameter space, these approximation methods are surprisingly 
accurate and can provide insight into the complex nature of the reversal process in ECC 
media. 
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A. Direct Path Approximation 


Figures [3] and [4] show that the energy calculated along the paths used as an initial guess in 
the determination of the MEP are in fact very close to those given by the MEP for both the 
strong coupling case (/ = 2.0 x 10 ~ 3 J/m 2 ) and the weak coupling case (/ = 0.5 x 1CP 3 J/m 2 ). 
This suggests that, in the strong coupling case, a reasonably good approximation to the rate 
coefficients can be found by replacing the MEP with the direct path 0\ = 6^ = 6. For this 
path the energy can be written as, 

E = - (K a v a + K b v b ) sin 2 6 - fi 0 H (M a v a + M b v b ) sin 6 - I A. (15) 


This expression for the energy is of the SW form and hence the expressions for the attempt 
frequency and energy barrier can be found analytically using the expressions in Brown’s 
classic paper , 2 


fa/3 ~ 


KtV i tto 7 
7ik B T v 1 + Oil 


H 2 

H k 2 


(H k ± H ), 


AE a s = —K t v I 1 =F —— 


H 


K 


(16) 

(17) 


where H K = 2 K T /M T , K T = K a + K b , and M T = M a + M b . In Fig. [9ja), we show a 
comparison of the MH hysteresis loops calculated using the rate coefficients calculated using 
the MEP method and the direct path approximation, with «o = 0.1 and J=2.0xl0 ~ 3 J/m 2 
for several sweep rates. 

Similarly in the weak coupling case, we can replace the four MEPs that link the minima 
{(Ti tA 02,02 cr 3 ,o - 3 tA cr 4 , } by the paths {6 ai 6 b j e {{ 6 », —tt/ 2 }, {tt/ 2 , 6 >}, {-vr/ 2 , 6 »}, { 6 », tt/ 2 }}. 

ft is straightforward to show that along each of the paths, the energy will be of the SW 
form and the rates may be calculated using Eqs. (16) and ( |I7| ). In Fig. |9](b), we show 
a comparison of the MH hysteresis loops calculated by the MEP method and the direct 
path approximation with J=0.5xl0 ~ 3 J/m 2 , for several sweep rates. The coercive field is 


shown as a function of sweep rate for 1= 2.0, 0.5, and O.lxlO -3 J/m 2 in Fig. 9(c), using both 
methods. As can be seen, for both the strong and the weak coupling cases, the differences 
between the MH hysteresis loops calculated from the MEP (exact) formulation and direct 
path approximation are generally small and only weakly dependent on the sweep rate R. 

One drawback of this approach is the fact that it is actually two distinct approximations, 
one valid for the strong coupling regime and another valid for the weak coupling regime, 
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FIG. 9. (a) shows the MH hysteresis loops at different sweep rates for 7=2.0x10 3 J/m 2 and (b) 
for 7=0.5 xlO 3 J/m 2 . Solid lines are obtained from the MEP method and the dashed lines are 
from the direct path approximation. The extracted coercivity as a function of the sweep rate for 
7=2.0, 0.5, and O.lxlO -3 J/m 2 is shown in (c). 


and it does not really provide an obvious way of interpolating between them. 


B. Transient State Approximations and Metabasins 


The second approximation to consider is based on the fact that, depending on the param¬ 
eters, there can be significant differences in the energy barriers and the attempt frequencies 
separating the energy minima. By way of an example, the calculated energy barriers and 
attempt frequencies between minima are presented in Table [I] together with the calculated 
rate constants r a p and the mean escape times r a p = 1 /r a p for the case 7 = 0.5 x 10 ~ 3 J/m 2 


and H = 0, shown in Fig. 2(a) Because of the factor exp(— AE/ksT) in the Arrhenius- 


Neel expression, the differences in A E a p (which are approximately 4 ~ 5) can lead to rate 
coefficients that differ by several orders of magnitude. This suggests that some of the states, 
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AEup/ksT 

f a p (GHz) 

r af 3 (MHz) 

T a p (PS) 

1 2 

12.4826 

20.5832 

7.80502 x 10“ 2 

1.28121 x 10 1 

2 -> 1 

3.79107 

8.47542 

1.91303 x 10 2 

5.22731 x 10“ 3 

1 -> 3 

19.6502 

51.7315 

1.51275 x 10“ 4 

6.56705 x 10 3 

3 ->• 1 

10.9587 

21.3012 

3.7078 X 10" 1 

2.69701 

2 —>■ 4 

10.9846 

8.52533 

3.60141 x 10" 1 

2.77668 

4—^2 

19.6762 

51.5684 

1.46935 x 10 -4 

6.75972 x 10 3 

3 —>■ 4 

3.78935 

8.52533 

1.92759 x 10 2 

5.18781 x 10 -3 

4^3 

12.4809 

20.7044 

7.86445 x 10“ 2 

1.27153 x 10 1 


TABLE I. The energy barriers, attempt frequencies, rate coefficients and mean escape times cal¬ 
culated from the MEPs connecting the minimum energy states for the case I = 0.5 x 10~ 3 J/m 2 


and H = 0 corresponding to the energy landscape shown in Fig. 2(a) 


in this case specifically states oy and cr 3 , are very short lived and will not contribute sig- 
nihcantly to the magnetization for processes involving long time scales (i.e. MH hysteresis 
loops generated using the vibrating sample magnetometer (VSM)). However, one has to be 
careful in removing such transients as they serve as intermediate states in the process of 
grain reversal. 

Consider for example state ay with both grains aligned along the positive z-axis. It can 
make the transition to states oy or oy. Comparing the mean escape times associated with 
the two transitions it is obvious, since 7 y _>. 2 <C Ti ^ 3 , that the predominant transition will be 
to state oy. From state oy the grain again has two choices. It can make the transition to 
state oy or back to state oy. Comparing the mean escape times it is clear, since t 2 ^i -C t 2 _> 4 , 
that the predominant transition is for the grain to return to its initial state ay. This implies 
that grains in the states o\ and oy will fluctuate back and forth with a characteristic time 
scale of the order of 10 /is for some time before it will transition to 1 —> 3 or 2 —» 4. The 
effect of these fluctuations will be to establish a local thermodynamic equilibrium between 
the two states oy and oy with a time scale on the order of a fraction of a ms. When local 
equilibrium is established, the net average probability flux between the two states will be 
zero and hence Ii _>. 2 = I 2 _>.i. A similar argument may be applied to the states oy and oy. 

The above argument implies that, while P\(t) and p 2 (t) are time dependant, they will 
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nevertheless satisfy the constraint, 


Pi(t) _ r 21 _ ( AG 12 \ 

P 2 (t) ~ r 12 ~ CXP { k B T ) 

where AGi 2 = G\ — G 2 and G a is expressed in terms of Z a , dehnecl in Eq. Q, as 


(18) 


G a = -k B T\ogZ a . 


(19) 


This is consistent with the requirement that states in the metabasin = fliUfl 2 satisfy the 
condition of local equilibrium Ci(t) = c 2 (t) = Ca (t) and hence the probability density within 
the metabasin formed by the union 72^ = U fl 2 is given by a Boltzmann distribution 
Pa{x, t) = c^(t) exp(— E(x, t)/k B T). Again a similar argument can be made for grains in 
the states er 3 and oq. 

The above reasoning implies that for processes with time scales on the order of ms or 
greater, we can assume that Pi(t)/p 2 (t) = r 2 i/ri 2 and p±(t)/pz(t) = r^A/r^. If we therefore 
define metabasins as those regions of phase space = hA U fl 2 and fl B = 0 3 U H 4 , then 
the probability of finding a grain in one of these metabasins is simply given by PA^t) — 
Pi(t) + p 2 (f) and p B {t) = p 3 (t) + p^ik) which can be shown to satisfy the following rate 
equations, 


dp A 

dt 

dp B 

dt 


~rABPA + r BA PB , 
~ r baPb + t ABPAi 


( 20 ) 


where the rate coefficients bab and r B A are given by, 


tab = 


r i 3 r 2 i + ri 2 r 24 


rBA = 


n 3 ir 43 + r 34 r 42 


ri2 + r 2 1 r 34 + r 43 

Using this concept of metastates, the set of four rate equations has been reduced to two, 
where Pa(H) and p B (H) can be obtained by numerical integration. The probabilities p a {H ) 
for a G {1,2, 3,4} can be determined from the values of Pa(H) and p B (H) together with 
the ratios r 21 (H)/r 12 (H) and r 43 (hf)/r 3 4 (ih) and hence the magnetization calculated as a 
function of H. 


Fig. 10 shows a comparison of the MH hysteresis loops for 7=0.5 x 10 ^ 3 J/m 2 between the 
original four-state model (solid curves) and the two-state approximation (open circles). The 
two models show very good agreement up to R < 10 s Oe/s above which the assumption of a 
Boltzmann probability distribution within the metastates A and B is no longer accurate. The 
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FIG. 10. A comparison of the MH hysteresis loops for 7=0.5x10 3 J/m 2 for the four-state model 
(solid lines) and the two-state approximation (open circles). 

above analysis in terms of metastates not only simplifies the system of equations that need 
to be solved for a range of R values but also provides some insight into how to understand 
the complex relationship between the sweep rate R and the response of ECC media. It is 
also important to note that while the case in which the system is described in terms of two 
metabasins, how the phase space up is divided into metastates for a given set of parameters 
is dependent on the nature of the energy landscape and time scales of interest. Indeed it is 
possible to adjust the number and regions of phase space occupied by the metabasins as the 
system evolves. In contrast to the previous approximation schemes, note that the two state 
model described in this section evolves smoothly into the coherent rotation of the strong 
coupling case as the exchange coupling constant / increases. 


VII. METASTATES AND KINETIC MONTE CARLO 

When the present model is extended to include magnetostatic and intralayer exchange 
interactions, the direct integration of the rate equations is no longer feasible. An alterna¬ 
tive approach is the KMC algorithm, which utilizes a stochastic algorithm to integrate the 
rate equations, and which can be adapted to include the interactions between the grains. 25 
However, the presence of low energy barriers can significantly increase the simulation time, 
effectively rendering the KMC approach no longer feasible at low sweep rates. This is 
a longstanding problem with KMC simulations.® 1 One way of dealing with this problem 
is by combining clusters of minimum energy states separated by low energy barriers into 
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“metabasins” as described in the previous section. 

To demonstrate the significance of the role of metabasins in the application of the KMC 
algorithm, consider the decay of an initially fully polarized ensemble of N identical non¬ 
interacting grains (pi = 1) with zero held and / = 0.5 x 10~ 3 J/m 2 . Using the rate coefficients 
presented in Table [T| the wait times for each of the N grains is given by 

t a ^ p (n) = r~l log(x), (21) 

where x is a uniformly distributed random number G{0<x<l},a denotes the state of 
the n th grain and f3 represents the two possible states it can transition into. The wait times 
describe how long we might expect to wait before the n th grain would make the transition 
a —> f3. The shortest of these wait times defines the first reversal time. The KMC step then 
takes the transition with the shortest reversal time and switches the grain from state a to 
a new state f3. This process is then repeated generating a stochastic sequence that models 
the process of thermally activated grain reversal. 



FIG. 11. Decay of the normalized magnetization at zero applied field with J=0.5xl(U 3 J/m 2 
(a) for the two-state model (red line) and the four-state model (blue line) from the KMC method 
together with the numerical integration of the four state model (black line) and (b) comparison of 
KMC results for two state model averaged over 1000 runs (red line) together with results obtained 
from numerical integration of the four state model (black line). 


Figure 11 (a) shows the magnetization m plotted as a function of t calculated using the 
KMC method for both the four state model and the two “metastate” representation for 
a system of 1000 non-interacting grains, together with the numerical solution of the rate 
equations for the four state model. The KMC solutions show the effects of the stochastic 
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fluctuations and both are in reasonable agreement with the solution obtained by direct 
integration of the rate equations. However, for the four state model, the average first reversal 
time was 6.1038 x lCU 9 s while for the two metastate representation, the average first reversal 
time was 3.1595 x 10~ 6 s; a factor of approximately 500 times greater than the four state 
case. This difference arises from the fact that vast majority of KMC steps in four state 
model were simply fluctuations within the metabasins A (ay <H> <72) and B (<73 <74). The 

difference in the average KMC time step is reflected in the run times; 39 minutes in the 
case of the four state model and approximately 4 seconds in the case of the two metastate 
representation. The speed up factor of 600 in completion times for the four and two state 
representations includes not only the shorter time steps but also the computational overhead 
associated with the four state model. To demonstrate the equivalence of the results obtained 
from the two state KMC calculations and those obtained by the direct integration of the 
four state model, Fig. [TT] (b) shows good agreement between a plot of the average mvs.t 
obtained from the two metastate representation averaged over 1000 independent KMC runs 
and those obtained by direct integration of the four state model. 

These results illustrate that for future applications with interacting grains, where direct 
integration of the rate equations is not feasible and the time scales of stochastic LLG re¬ 
stricts its application to p,s time scales, the KMC approach represents a viable model of 
long-time processes dominated by thermally activated reversal. Further, when the system 
in question, such as ECC media, has a range of energy barriers, the above example demon¬ 
strates that removing the short time scale fluctuations associated with transient states by 
combining them into a single metabasin can result in significant computational efficiencies 
with negligible loss of accuracy. In the case of interacting systems, the gains in run time 
are even more significant given the increased computational overhead involved in computing 
the effective fields due to the interactions and the more complex energy landscapes that 
typically include a greater number of critical points than the simple model discussed here. 


VIII. DISCUSSION AND CONCLUSIONS 

A set of rate equations are presented that describe the evolution of a non-interacting 
ensemble of dual layer ECC grains based on processes of thermally activated grain reversal. 
The rate coefficients are calculated from the Langer formalism and have the Arrhenius-Neel 


25 


form in which the attempt frequency and energy barriers are expressed in terms of the 
energy and its Hessian matrix calculated at the maximum point on the minimum energy 
paths that connect the energy minima. The particular form for the attempt frequency is 
outlined in the Appendix and is not restricted to the canonical coordinates commonly used 
in the derivation but is valid for any system (or systems) of generalized coordinates that 
parameterize the surface of a sphere. The minimum energy paths are calculated using the 
so called “string method”. The rate equations can be integrated numerically for the case of 
a time dependent applied held with a constant sweep rate and the magnetization calculated 
to produce MH hysteresis loops. 

It is shown that the method may be used to study both the strong coupling regime, 
in which the energy landscape has two energy minima, consisting of two ferromagnetically 
aligned layers as well as the more complicated weak coupling regime, which has an energy 
landscape that can have up to four distinct energy minima, two ferromagnetic and two 
antiferromagnetic states. Calculating the MH hysteresis loops therefore requires solving 
two coupled rate equations for the strong coupling regime and up to four coupled rate 
equations for the weak coupling regime. The results show that, for the parameters used in 
the current work, the transition from the weak to the strong coupling regime occurs when 
/ ~1.0—1.5xl0 _3 J/m 2 , which is the region of interest for ECC based recording media. 

Verification of our model results for MH hysteresis loop was achieved through comparison 
with LLG simulations on a dual layer system, each layer with a 16 x 16 non-interacting 
grains. The high degree of agreement confirms the accuracy of the rate coefficients and the 
numerical integration of the rate equations. In addition, using the MH hysteresis loops 
calculated from the rate equations, the effect of rate dependence and exchange coupling on 
a Figure of Merit based on the ratio between the energy barrier and switching held was 
calculated. This provides some guidance on the optimal coupling between the layers. 

Results based on a direct path approximation to the MEP in the strong and weak cou¬ 
pling limits that permit analytic expressions for both the energy barriers and the attempt 
frequencies are presented in Figs. |9|a) and[9](b). The results show remarkably good agree¬ 
ment with the exact MEP calculation for both the strong, / = 2.0 x 10 -3 J/m 2 , and the 
weak, / = 0.5 x 10 -3 J/m 2 , coupling regimes. 

Another approximation scheme was presented in which pairs of minima separated by a 
relatively low energy barrier so that they are very rapidly equilibrated and could be combined 
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into a single metastate in which the ratio p a /pp is given by a Boltzmann factor (Eq. (18)). 
It was shown that for the case / = 0.5 x 10 -3 J/m 2 , the MH hysteresis loops obtained 
by integrating the four state rate equations could be accurately reproduced by integrating 
the rate equations for a two “metastate” representation with rate coefficients between the 


metastates given by Eq. (20). The potential importance of this mapping of “exact” the 
four state model to a model consisting of two metastates was demonstrated in simulation of 
magnetic decay using the KMC algorithm, in which the two “metastate” model produced 
results essentially equivalent to the four state model, but with a run time that was reduced 
by a factor of 600. 

The results of this work serve as a prelude to the extension of our previous KMC 
approach^ 7 to study thermally activated magnetic grain reversal in dual layer ECC media 
that includes magnetostatic and intralayer exchange interactions) 25 The essentially exact 
treatment of grain reversal for the dual layer ECC grain problem as outlined in this work, 
serves as the foundation for this extension, while combining cluster of states that are sep¬ 
arated by relatively small energy barriers to form metastates allows us to deal with the 
phenomena of “stagnation” that can severely limit the accessible run times that can be 
achieved using the KMC approach. 

This extension of our previous KMC algorithm will allow for the direct comparison of 
experimentally determined slow-sweep-rate MH hysteresis loops for ECC media with corre¬ 
sponding modelled results. This capability is useful for the estimation of model parameters 
which characterize recording media such as intralayer and interlayer exchange couplings. 
Such a direct comparison is not possible with traditional LLG simulations where long time 
scales are inaccessible. This dual layer KMC algorithm will also be especially useful in 
applications to dual layer media for heat assisted magnetic recording where thermally acti¬ 
vated moment reversal is pronounced. 20 In addition, the investigation of magnetostatic and 
intralayer interaction effects on the FOM of Fig. [8] is of particular interest. 
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Appendix A 


The attempt frequency given in Eq. ([7]) is key to the analysis presented in the previous 
sections, we therefore outline the derivation in some detail. The approach starts with the 
Fokker-Planck equation (FPE) for a single magnetic moment in a magnetic field and the 
generalization to consider a set of exchanged coupled moments. 

The FPE can be derived from the stochastic LLG equation. 22 29 30 From the FPE, the rate 
constants r a p defined by Eq. ([6]) are calculated using the formalism presented by Langei 22 
adapted to account for the dissipative dynamics of magnetic moment in an applied field. 32 
While the application of the Langer formalism is facilitated by a judicious choice of coordi¬ 
nates, (0, z = cos 6), often referred to as the canonical coordinates, it is nevertheless possible 
derive a straightforward expression for the rate coefficients based on any set of generalized 
coordinates (w 1 , u 2 ) that parameterize the surface of the sphere § using as basis vectors the 
covariant tangent vectors gj = drh/du l . An advantage of this approach is that it allows the 
direct application of the tools of differential geometry to be applied to the problem. This 
is of some practical importance in the case of spin dynamics as it is not possible to define 
a single coordinate system on the surface of a sphere where the metric is everywhere finite. 
However, the surface of a sphere can be treated as a differentiable manifold by dividing it 
into overlapping regions, each of which has a metric that is everywhere finite. 

Consider a magnetic moment m a of volume v a with magnetization M a , anisotropy con¬ 
stant K a , and a damping factor «o in a magnetic field H. The equation of motion for the 
moment is given by 


drh a 

dt 


Wo m a 


x H 


Oi 0 


1 + Q! 2 


0 


rh n x m n x H 


(Al) 


where rh a = m a /M a v a and H = —Hq 1 dE/drh a . We parameterize the unit vector rh a in 
terms of the generalized coordinates u = (m 1 ,^ 2 ) (e.g. u = (0 a , </>„)) which cover the surface 
of the unit sphere. Since drh a /dt will be tangential to the surface of the sphere § a , we define 
the local covariant basis vectors® 


_ dm a 
^ did 
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(A2) 







Any vector tangential to the surface of the sphere can therefore be written as v = giV 1 , where 
the components v l define a type (1,0) tensor. The basis vectors < 7 * are, in general, neither 
orthogonal nor normalized to unity, but satisfy 

9i ' 9j = 9iji (A3) 

where g tl is the metric tensor. We also define the reciprocal, or contravariant, basis vectors 
<f such that 

Si ■ S’ = H, (A4) 

where 5j is the Kronecker delta function. Any tangential vector v may then also be written 
in contravariant form as 


v = g Vi. 

The scalar product of any two tangential vectors u and v may then be written as 

u ■ v — Ui.v 1 = g^UiVj = g.ijU l v\ 


(AS) 


(A 6 ) 


where the components v t define a type (0,1) tensor. The vector rh a may also be written in 
terms of the basis vectors g t and g l as 

9i x 92 JiX g 2 


rrin = 


. 9*1 x g 2 \ | g 1 x g 2 \' 


It is straightforward to show that 


rh a x gi = \fg€ij g J 




m a x g = ~j=9v 

where e V] and e 13 denote the Levi-Cevita symbols defined as 


e n — ^22 — 0, 612 — —£21 — 1, e 11 — — 0, e 1 " — —e zl — 1, 


,22 


,12 


,21 


(A7) 

(AS) 

(A9) 

(A10) 


and g = m a ■ (g 1 x g 2 ), which is simply the volume of the vectors triad (m 0 , gi, g 2 ). It can 
be shown that e 13 /^/Jj and \fge^ define tensors of the form ( 2 , 0 ) and ( 0 , 2 ), respectively, on 
the surface §. The LLG equation can then be written in covariant form as 


du l 


( e ij 


v = ~nr = -7s - 


dt 

JB_ 

m a \y/g 1 + a§ : 






Wg 1 + «r J 

dE(u ) 


1 do Hj 




duj 


(All) 
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We note that the form of the above equations are invariant under a generalized coordinate 
transformation. 

Consider an ensemble of such spins and denote by p{u, t ) the probability density, then 
the probability of a single spin will be aligned in the solid angle dVt at time t is given by 
dp(u,t) = p{u,t)dVL. The probability density p(u,t) may be calculated from the Fokker- 
Planck equation (FPE) which can be written in terms of the coordinates u l as 


dp(u, t) 
dt 


-Vij{u,t), 


(A12) 


where the probability current density J* (u, t) consists of an advective term and a diffusive 
term 


J‘ (u, t) = p(u, t)v‘ - dsBlLgV V . p („, t), 

i a 0 

where D a = aok B T/^ B m a , with m a = M a v a , the velocity held v l is given by Eq. 
V, denotes the absolute derivative, and 


(A13) 


(All) and 


Vipiu) 

Viv\u) 


dp(u ) 
dui 

1 d(v^t/(n)) 
dui 


dv l (u) 

du l 


+ v i (u) V, 


(A14) 

(A15) 


where Y t j k denotes the Christoffel symbol of the second kind. 33 Since we are interested in 
solutions close to equilibrium, following Langerj- 22 ^ we write the probability density in terms 
of the crossover function c(u, t ) as 


p(u,t) 


c(u, t ) exp 


f Eg(U,t) \ 

V k B T )- 


(A16) 


It can then be shown that J L (u,t) may be written in terms of the crossover function c(u,t) 
as 


J l (u,t) =&BT^-ex p 
m„ 


E a (u, t) 




«o 


k B T ) \y/g 
divergenceless terms. 


i 9«- 

l-a 2 0 


g lJ Vjc(u, t ) 


(Air) 


The above formalism can be readily extended to the problem of two coupled spins. Let 
w = (w 1 ,™ 2 ) denote the generalized coordinates that specify the orientation of a second 
magnetic moment of volume Vb, magnetization Mb, anisotropy constant K h and damping 
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constant ao- In the absence of the interaction, the probability current density on the surface 
of the sphere S 5 may then be written as 


T if ( E b (w,t)\ (e 3 a 0 , , 

J w,t =k B T —exp- —— -= - ——j 9 Vjc(w,t) 

m b V k B T J Vvy 1 + ag ' 

+ divergenceless terms. 


(Ais) 


In presence of an interaction between the moments, we define the vectors x = (u, w ) that 
spans the tangent space of the four dimensional manifold § = § a <E> §&. The metric g^ v can 
be written in matrix form as 

V'MII 0 

0 Il9 y (i>)ll 

where ||(f J (a)|| and \g l3 (b) \ | denote the matrix forms for the metrics in the manifolds § a and 
§& for the single grains a and b. The LLG equation for the case of interacting spins may 
then be written in covariant form as 


ll<rii = 


(A19) 


dx 11 = 7 B dE(x,t) 

dt m dx v 

where m = m a + rn b and the tensor T tiU (x) expressed in matrix form as 


(A 20 ) 


/ 


117 ^ 11 — 


m 

m n 


Oi 0 


1 + Q?n 


llr'WIH 


f v n 


\fK^ 


m 


ot 0 


||e^|| \ 


V 


m„ \ 1 + 4 I|93(6)I1 yW)JJ 


(A21) 


This yields the following expression for the probability current density in terms of the 
crossover function c(x, t) 

= - k B T 0) exp T^(x)V„c(x,t) 

+ divergenceless terms. (A 22 ) 

This gives 


V^J M (a;,f) = -k B T—ex p f- 
m \ 


E 


«o „ „ uu 1 dE 


k B T J Vl + « 2 M 




k B T dx., 


T^)\7 u c{x,t), (A23) 


where G^ may be written in matrix form as 



(— irwii 

0 

\ 

II = 

m a 

— II 9 ij 



l 0 

(&)ll, 


V 

m h 


(A24) 
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As discussed in Secs. IIA and IIB we are interested in stationary solutions that satisfy 


= 0 for which the crossover function is essentially homogeneous except in a narrow 
region in the neighbourhood of the boundaries r a/ 3 where it goes from c Q —> cp on crossing 
the boundary from D Q —>• VLp. These solutions correspond to a state of “local” equilibrium 
with thermodynamic equilibrium corresponding to the special case c a = const for all a. In 
addition, as discussed in Sec. |TT| for the energy scales we are interested in, the probability 
current density is concentrated in a narrow region surrounding the saddle point s a p on the 
boundary Y a p. The crossover function is thus required only in region surrounding s a p. This 


allows for two approximations that simplify Eq. (A22). The first is to assume that the 


coordinate system is chosen such that the metric does not have any singularities close to 
the saddle point and it can be approximated as a constant. The second, assumes a quadratic 
approximation for the energy 

d 2 E(x) 


E(x) « E(x s ) + - ^ 


(x - x s y{x - x s ) v + 


2 ^ dx^dx” 

Defining the eigenvectors and eigenvalues of the Hessian matrix d 2 E(x) / dx^ L dx u \ x=Xa as 

d 2 E(x 


(A25) 


dx^dx v 


< = A n <, 


(A26) 


we define the new coordinates y n as 


y n = a l ( x ~ x sT , 


(A27) 


with a™a v n = S rnn . The quadratic form of the energy may then be written as 


E(x)^E(x,) + -Y,^n(y n ) 2 + --- 


(A28) 


71=1 


Note that Ai > A 2 > 0, A 3 = 0 (by symmetry), and A 4 < 0. From Eq. (A24) we then obtain 
in the static limit (V M J M (x) = 0 ) the following equation for the crossover function 


V' (A-G™ _ J_ 8 -M = 0, 

^ U +aldy m k B T u 1 J -" 


m,n 


dy r 


(A29) 


where G mn = T mn = and omits the term rn — 3 and n = 3. As 

discussed by Langer, 22j this equation may be solved using the method of characteristics, 
whereby we look for solutions of the form c(y) = c(t ) where the variable t defines a trajectory 
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t = Yhn .^3 U n y n , with the direction cosines U n are given by the solutions of the eigenvalue 
equation 


X m T mn U n = £U m . (A30) 


The solution of interest is given by 


dc(t) 
dt 



(A31) 


with k = £(1 + aD/aoGuksT, where £ denotes the negative eigenvalue obtained form 


Eq. (A30) with Gjj = G mn U m U n . Integrating this equation using the boundary conditions 
lim^.oo c(t) = p a /Z a and lim^oo c(t) = Pp/Zp gives 


dc(t ) 
dt 

Writing t in terms of the direction cosines U n gives dc(y)/dy n = U n dc/dt leads to the 
following expression for the probability current density in the region around the critical 
point 


'l«l lPy_ 

27T \Zp 


Pa 

Z n 


exp 


\hi\t 2 


(A32) 


J m (y ) 


«o 


1 + CTq 

Urr 

Xrr 


G 


u 




, n rp\2 f Pp Pa\ ( E s 

^ (ksT) \Tfi ~ ~zj exp \kPT 


\ nk 


Az- 


k B T 


&nk A K | U n U k ) PnVk 


(A33) 


for m pX 3 (J 3 = 0). 

To calculate the net probability current Z a ^.p flowing between the basins of attractions 
and flp, we simply integrate the T 4 component of probability current density over the 
hypersurface defined by t/ 4 = 0 to give 34 




VgJy)J\y) 


y 4 =o 


dy 1 dy 2 dy 3 , 


(A34) 


where g s (y ) is dehned in terms of the metric associated with the subspace formed by the 
vectors {y 1 , y 2 , y 3 } 


9u(y) 9i2(y) gi3(y) 


9 s(y) = det 


hi{y) 922(y) 923 ( 9 ) 


9zi(y) 93-2(9 ) 933(9) 


(A35) 
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Because of the exponential factor in the expression for the probability current density, 


Eq. (A33), only the region in the immediate vicinity of y 1 = y 2 = 0 will contribute to the 


integral and we can therefore use the quadratic form of the energy given by Eq. (A25). Also 
it is convenient to choose y 3 so that it corresponds to the azimuthal angle <f> = (4> a + 4>b )/2 as 
the integration with respect to 7/3 simply yields a factor of 2n. The net probability current 
l a ^p may then be evaluated to give 


i a ^ = VW)-g 


7 B r, « 0 

U 

m 


1 + (Un 




(2vr k B Tf 


AoA 


exp 


2 / ' 41 


E x 


k B T 


Writing the net probability current as I a ++p = — Z a ^p yields 


I a ^ = -\f?j{s)—G 
m 


1 + OL Q 


(27 Tk B Tf 

IA1A2A4 


exp 


E„ 


k B T 


(A36) 


(A37) 


and hence the following expression for the rate constants r a p 


r a/3 


= VW)~ G 


«o 


UT~. -2 

m 1 + ttn 


\K\ 


(exp (- E s /k B T) \ / (2irk B T) 3 




Z n 


IA1A2A4 


(A38) 


In order to compute Z a , again assume that the probability density is strongly localized 
at cr a and, since Z Q is a scalar quantity, transform from the coordinates x M to some new 
coordinates x ^ so that the metric g^ has no zeros or singularities in the region of interest. 
Thus 

'd 2 E(x)' 


Z a = a/ g(oi)(2nk B T) 2 det 


I g(a)(2nk B T) /l 




dx^dx v 

exp 


exp 


E a 

ErT 


x=x (a) 


E n 


k B T 


where ry denotes the eigenvalues of the Hessian matrix d 2 E(x)/dx^dx v \- = 
Eq. (A40) into Eq. (A38) gives the result in Eq. Q 

71727374 _ ( (E s - E, 


x=x(a)' 


T aft = 


«o g(s ) 7b 

1 + c*o V di a ) m 


Gu\k\ 


2tt k B T | A 1 A 2 A 4 1 


exp 


kuT 


(A39) 

(A40) 

Substituting 

(A41) 
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( 2011 ). 
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